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Abstract 

In the past few years a lot of attention has been paid in the multigrid 
solution of multilevel structured (Toeplitz, circulants, Hartley, sine (r 
class) and cosine algebras) linear systems, in which the coefficient matrix 
is banded in a multilevel sense and Hermitian positive definite. In the 
present paper we provide some theoretical results on the optimality of 
an existing multigrid procedure, when applied to a properly related al- 
gebraic problem. In particular, we propose a modification of previously 
devised multigrid procedures in order to handle Hermitian positive def- 
inite structured-plus-banded uniformly bounded linear systems, arising 
when an indefinite, and not necessarily structured, banded part is added 
to the original coefficient matrix. In this context we prove the Two-Grid 
method optimality. 

In such a way, several linear systems arising from the approximation of 
integro-differential equations with various boundary conditions can be ef- 
ficiently solved in linear time (with respect to the size of the algebraic 
problem). Some numerical experiments are presented and discussed, both 
with respect to Two-Grid and multigrid procedures. 

1 Introduction 

In the past twenty years, an extensive literature has treated the numerical so- 
lution of structured linear systems of large dimensions [11) . by means of pre- 
conditioned iterative solvers. However, as well known in the multilevel setting, 
the most popular matrix algebra preconditioners cannot work in general (see 
P51 EH |21] and references therein), while the multilevel structures often are 
the most relevant in practical applications. Therefore, quite recently, more at- 
tention has been focused (see [I] [21 HH1 HS1 H31 HE1 HH1 EZ1 W\) on the multigrid 
solution of multilevel structured (Toeplitz, circulants, Hartley, sine (r class) and 
cosine algebras) linear systems, in which the coefficient matrix is banded in a 
multilevel sense and Hermitian positive definite. The reason is due to the fact 
that these techniques are very efficient, the total cost for reaching the solution 
within a preassigned accuracy being linear as the dimensions of the involved 
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linear systems. 

In the present note we propose a slight modification of these numerical multigrid 
procedures in order to handle structured-plus-banded uniformly bounded Her- 
mitian positive definite linear systems, where the banded part which is added 
to the structured coefficient matrix is indefinite and not necessarily structured. 
A theoretical analysis of the related Two-Grid Method (TGM) is given in terms 
of the algebraic multigrid theory considered by Ruge and Stiiben [25]. More 
precisely, we prove that the proposed TGM is optimally convergent with a 
convergence rate independent of the dimension for a given sequence of linear 
systems {B n y n = c„}„ with uniformly bounded Hermitian positive definite 
matrix sequence {B n } n , under the assumption that such TGM is optimal for 
n%n — ^n}n with a given Hermitian positive definite matrix sequence {^4n}n 
related to {B n } n by means of a simple order relation. More precisely, we require 
A n < $B n , with || -Bn|| 2 < M, for some i9, M > independent of n and for every 
n large enough. 

As a case study, we may consider the case where B n = A n + 0„ where A n is 
structured, positive definite, ill-conditioned, and for which an optimal multigrid 
algorithm is already available, and where O n is an indefinite band correction 
not necessarily structured; moreover, we require that {^4 n } n and {B n } n are 
uniformly bounded and that A n + 0„ is still positive definite and larger than 
A n /d for some d > independent of n. 

For instance such a situation is encountered when dealing with standard Finite 
Difference (FD) discretizations of the problem 

Cu = — V 2 u(a;) + [i{x)u{x) — h(x), x G fi, 

where n{x) and h(x) are given bounded functions, Q, = (0, l) d , d > 1, and with 
Dirichlet, periodic or reflective boundary conditions (for a discussion on vari- 
ous boundary conditions see [HJ 29 ). For specific contexts where structured- 
plus-diagonal problems arise refer to [HI HHj and [5] , when considering also a 
convection term in the above equation. However, the latter is just an example 
chosen for the relevance in applications, but the effective range of applicability 
of our proposal is indeed much wider. 

The numerical experimentation suggests that an optimal convergence rate should 
hold for the MGM as well. Here, for MGM algorithm, we mean the simplest 
(and less expensive) version of the large family of multigrid methods, i.e., the 
V-cycle procedure: for a brief description of the TGM and of the V-cycle algo- 
rithms we refer to Section [H while an extensive treatment can be found in [17] 
and especially in [34] . Indeed, we remark that in all the considered examples 
the MGM is optimal in the sense that (see [1]): 

a. the observed number of iterations is constant with respect to the size of the 

algebraic problem; 

b. the cost per iteration (in terms of arithmetic operations) is just linear as 

the size of the algebraic problem. 

Nevertheless, it is worth stressing that the theoretical extension of the optimality 
result to the MGM is still an open question. 

The paper is organized as follows. In Section [2] we report the standard 
TGM and MGM algorithms and we write explicitly the related iteration matri- 
ces. In Section [3] we first recall some classical results related to the algebraic 
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TGM convergence analysis and then we prove the optimal rate of convergence 
of the proposed TGM, under some general and weak assumptions; the MGM 
case is briefly discussed at the end of the section. In Section [3] we analyze in 
detail the case of the discrete Laplacian-plus-diagonal systems and in Section 
we report several numerical experiments. Lastly, Section [6] deals with further 
considerations concerning future work and perspectives. 

2 Two-grid and Multigrid Method 

Let no be a positive d-index, d > 1, and let N(-) be an increasing function with 
respect to uq. In devising a TGM and a MGM for the linear system 

^no^no — ^noj (2T) 

where A no G C Ar( "« )x,v(ni,) and x no ,b no G C N ( n °\ the ingredients below must 
be considered. 

Let rii < uq (componentwise) and let € QN(n )xN{ni) ^ e a gj ven f u ll_r an k 
matrix. In order to simplify the notation, in the following we will refer to any 
multi-index n s by means of its subscript s, so that, e.g. A s :— A ria , b s := b ria , 
P s s +1 := Pnl +1 , etc. With these notations, a class of stationary iterative methods 
of the form 

x^+V = V s x{ j) + b s 

is also considered in such a way that Smooth(xi\ b s ,V s ,i> s ) denotes the appli- 
cation of this rule v s times, with v a positive integer number, at the dimension 
corresponding to the index s. 

Thus, the solution of the linear system (|2.ip is obtained by applying repeatedly 
the TGM iteration, where the j th iteration 

Xq +1 = TGM(Xq \ bo, Aq, Vb,pre) ^O.prej Vb.post) ^Ojpost) 

is defined by the following algorithm [T?] : 

yo ■= TGM(x ,b , A , V , pie , ^o.pre, Vq ,post ) ^0,post ) 



1. x ■= Smooth(x ,b , V 0tPie , v , pTe ) 

2. r := b a - A x 

3. n := {pl) H r 

4. Solve A x yi = r u with A 1 := {pl) H A pl 

5. y := io+PoVi 

6. y := Smooth(y ,bo, V ,post ) *^0,post ) 
Steps 1. and 6. concern the application of i/q pre 

steps of the pre- smoothing (or 
intermediate) iteration and of vo tPO st steps of the post-smoothing iteration, re- 
spectively. Moreover, steps 2. — > 5. define the so called coarse grid correction, 
that depends on the projection operator (pq) h ■ In such a way, the TGM itera- 
tion represents a classical stationary iterative method whose iteration matrix is 
given by 

TGMo = < P os°r CGCo . (2-2) 

where 

CGCo =h-pl [(pI) H AopI\ -1 [pI) h Ao 
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denotes the coarse grid correction iteration matrix. 

The names intermediate and smoothing iteration used above refer to the multi- 
iterative terminology [26) : we say that a method is multi-iterative if it is com- 
posed by at least two distinct iterations. The idea is that these basic components 
should have complementary spectral behaviors so that the whole procedure is 
quickly convergent. In our case the target of the smoothing iteration is to reduce 
the error in the subspace where A$ is well-conditioned, but such an iteration will 
be slowly convergent in the complementary space. The coarse grid correction it- 
eration matrix is a (non-Hermitian) projector (see e.g. |27j ) and therefore shows 
spectral radius equal to 1 . As a consequence, the corresponding iterative proce- 
dure does not converge at all, but it is very quickly convergent in the subspace 
where Aq is ill-conditioned, if is chosen in such a way that its columns span 
a subspace "close enough" to the ill-conditioned one. Finally, the intermediate 
iteration is strongly convergent in that subspace where the combined effect of 
the other two iterations resulted to be less effective. Notice that in our setting 
of Hermitian positive definite and uniformly bounded sequences, the subspace 
where Aq is ill-conditioned corresponds to the subspace in which Aq has small 
eigenvalues. 

Starting from the TGM, we introduce the MGM. Indeed, the main difference 
with respect to the TGM is as follows: instead of solving directly the linear sys- 
tem with coefficient matrix A\, we can apply recursively the projection strategy 
so obtaining a multigrid method. 

Let us use the Galerkin formulation and let uq > n\ > . . . > rt; > 0, with I being 
the maximal number of recursive calls and with N(n s ) being the corresponding 
matrix sizes. 

The corresponding multigrid method generates the approximate solution 

X^ + 1) = MGM(0,X%',b ,A , Vb,prc,fO,prc, Vb.post^Cpost) 

according to the following algorithm: 

Us := fiAGM[s, X S: b S: A s , V s ,p rCl ^s.prci Kg, post i ^s.post ) 

If s = I then Solve(A s y s = b s ) 

else 1. i s := Smooth (x s , b s , F s , pre , ^ s , pr o) 
.— h s - A s x s 



2. r s 

3. r s+ i 
4- Vs+i 



= (p s s +1 ) H r s 



:= MGM(s + l,0 s+ i,6 s+ i, A s+1 ,V s+ltPrc ,v s+ i, pic , 

Ks+l.post, ^s+l.post) 

5. y s := x s +pl +1 y s+ i 

6. y s := Smooth (y s ,b s ,V s , pos t,is StPOSt ) 

where the matrix A s+ i := (p s s +1 ) H A s p s s +1 is more profitably computed in the 
so called pre- computing phase. 

Since the multigrid is again a linear fixed-point method, we can express x^ 1 ^ 
as MGMqx^ +(7 - MGM )AQ 1 b , where the iteration matrix MGM is 
recursively defined according to the following rule (see [34]): 

(MGMi = O, 

MGM. = V%££ [l s -p s s +1 (I s+ i-MGM s+1 )A-i 1 (p s s +1 ) H As] V S %T, 

s = 0,...,l-l, 

(2.3) 
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and with MGM S and MGM s+ \ denoting the iteration matrices of the multigrid 
procedures at two subsequent levels, s = 0, ... ,1—1. At the last recursion level 
I , the linear system is solved by a direct method and hence it can be interpreted 
as an iterative method converging in a single step: this motivates the chosen 
initial condition M GMi = O. 

By comparing the TGM and MGM, we observe that the coarse grid correc- 
tion operator CGC S is replaced by an approximation, since the matrix A~^ 
is approximated by (I s +i — MGM s+ \) A~ +l as implicitly described in (|2.3[) for 
s = 0, 1. In this way step 4., at the highest level s = 0, represents an 
approximation of the exact solution of step 4. displayed in the TGM algorithm 
(for the matrix analog compare (|2.3p and lj2.2JI ). 

Finally, for I = 1 the MGM reduces to the TGM if Solve(A iyi = &i) is 
yi = Ai%. 

3 Discussion and extension of known conver- 
gence results 

Hereafter, by || • H2 we denote the Euclidean norm on C m and the associated 
induced matrix norm over C mxm . If X is Hermitian positive definite, then its 
square root obtained via the Schur decomposition is well defined and positive 
definite. As a consequence we can set || • \\x = WX 1 / 2 ■ \\2 the Euclidean norm 
weighted by X on C m , and the associated induced matrix norm. 
In the algebraic multigrid theory some relevant convergence results are due to 
Ruge and Stiiben [25]. In fact, they provide the main theoretical tools to which 
we refer in order to prove our subsequent convergence results. 
More precisely, by referring to the work of Ruge and Stiiben [25J , we will con- 
sider Theorem 5.2 therein in its original form and in the case where both pre- 
smoothing and post-smoothing iterations are performed. In the following all 
the constants a, a prc , a pos t, and f3 are required to be independent of the actual 
dimension in order to ensure a TGM convergence rate independent of the size 
of the algebraic problem. 

Theorem 3.1 Let Aq be a Hermitian positive definite matrix of size N(jiq), let 
p\ G (^N(n )xN(n 1 ) ^ no ^ n ^ ^ g fl gi ven full-rank matrix and let Vo. pos t be the 
•post- smoothing iteration matrix. 

Suppose that there exists a pos t > 0, independent of no, such that for all x £ 

\\V , poBt xf Ao < ||z||i -a post IMI^-i^, (3.1) 

where Dq is the diagonal matrix formed by the diagonal entries of Aq . 
Assume, also, that there exists (3 > 0, independent of uq, such that for all 

x € C N( - no) 

min \\x-ply\\l <(3 \\x\\\. (3.2) 

Then, [3 > a pos t md 

\\TGM q \\a < ^l-a post //3 < 1. (3.3) 
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Theorem 3.2 Let A$ be a Hermitian positive definite matrix of size N(no), let 
Pq € (£N(n )xN(n 1 ) ^ n ^ ^ n ^ fo e a gi ven full-rank matrix and let Vo, pre ; Vo jPOS t 
be the pre-smoothing and post- smoothing iteration matrices, respectively. 
Suppose that there exist a prc , a pos t > 0, independent of uq, such that for all 

x e c N ( n °) 

\\V 0A , rc x\\ 2 Ao < \\x\\ 2 A() - a prc \\Vo,pr C x\\ 2 AgD -i Ao , (3.4) 
ll^post^Hio < \\x\\ A() ~ a post \\x\\ 2 AoD - lAo , (3.5) 

where Dq is the diagonal matrix formed by the diagonal entries of Aq . 
Assume, also, that there exists (3 > 0, independent of uq, such that for all 

x e c N ( n °) 

\\CGC x\\ 2 Ao <f3 M\ 2 AqD -, Aq . (3.6) 

Then, (3 > a pos t and 

\\TGM \\ Ao < Jfc^£ < 1. (3-7) 
y 1 + a pro //5 

Remark 3.3 Theorems 13.11 and 13.21 still hold if the diagonal matrix Dq is re- 
placed by any Hermitian positive matrix Xq (see e.g. [2]). More precisely, 
Xq = I could be a proper choice for its simplicity, since any contribution due 
to the use of a different matrix will be subject to a formal simplification in the 
quotients a pre /[3 and ct pos t/f3. 

Remark 3.4 For reader convenience, the essential steps of the proof of Theo- 
rems 13.11 and 13.21 are reported in Appendix\^ where relations (|3.ip and (|3.4|) are 
called post- smoothing and pre-smoothing property, respectively, and the relation 
(|3.6p is called approximation property. In this respect, we notice that the ap- 
proximation property deduced by using (|3.2[) holds only for vectors belonging to 
the range of CGCq, see (|A.1[) : conversely the approximation property described 
in (|3.6[) is unconditional, i.e., it is satisfied for all x € C^*-™ -*. 

In this paper we are interested in the multigrid solution of special linear 
systems of the form 

B n x = b, B n e C N(n)xN(n \ x, b € C N(n) (3.8) 

with {B n } n Hermitian positive definite uniformly bounded matrix sequence, n 
being a positive d-index, d > 1 and N(-) an increasing function with respect 
to it. More precisely, we assume that there exists {A„}„ Hermitian positive 
definite matrix sequence such that some order relation is linking {A„}„ and 
{Bn}n, for n large enough and we suppose that an optimal algebraic multigrid 
method is available for the solution of the systems 

A n x = b, A n eC N ^ xN ^ n \ x,beC N( - n \ (3.9) 

We ask wether the algebraic TGM and MGM considered for the systems p.9[) 
are effective also for the systems (I3.8[) . i.e., when considering the very same 
projectors. Since it is well-known that a very crucial role in MGM is played by 
the choice of projector operator, the quoted choice will give rise to a relevant 
simplification. The results pertain to the convergence analysis of the TGM and 
MGM: we provide a positive answer for the TGM case and we only discuss the 
MGM case, which is substantially more involved. 
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3.1 TGM convergence and optimality: theoretical results 

In this section we give a theoretical analysis of the TGM in terms of the alge- 
braic multigrid theory due to Ruge and Stiiben |25j according to Theorem 13. II 
Hereafter, the notation X < Y, with X and Y Hermitian matrices, means that 
Y — X is nonnegative definite. In addition, {X n } n , with X n Hermitian positive 
definite matrices, is a uniformly bounded matrix sequence if there exists M > 
independent of n such that ||X„||2 < M, for n large enough. 

Proposition 3.5 Let {A n } n be a matrix sequence with A n Hermitian positive 
definite matrices and let Po € QN(n )xN(ni) ^ e a gj ven full-rank matrix for any 
uq > such that there exists Pa > independent of no so that for all x £ (C N ^ n °^ 

min \\x-ply\\ 2 2 <p A \\x\\ Ao . (3.10) 

Let {B n } n be another matrix sequence, with B n Hermitian positive definite ma- 
trices, such that A n < i3B n , for n large enough, with $ > absolute constant. 
Then, for all x £ <C N ^ n °^ and no large enough, it also holds [3b — I3a$ and 

min \\x - plyg < p B \\x\\ 2 Ba . (3.11) 

Proof. From (|3.10p and from the assumptions on the order relation, we deduce 
that for all x € C N{n) 

min \\x-plv\\l < 0a\\x\\ 2 Aq < ^ A \\xf Bo , 

i.e., taking into account Remark 13.31 the hypothesis (|3.2p of Theorem 13.11 is 
fulfilled for {B n } n too, with constant [3b — /3a$, by considering the very same 
projector considered for {A n } n . • 



Thus, the convergence result in Theorem 13.11 holds true also for the matrix 
sequence {£?„}„, if we are able to guarantee also the validity of condition (|3.ip . 
It is worth stressing that in the case of Richardson smoothers such topic is 
not related to any partial ordering relation connecting the Hermitian matrix 
sequences {v4„}„ and {B n } n . In other words, given a partial ordering between 
{A„} n and {B n } n , inequalities (|3.1|) . (|3.4p . and (|3.5p with {B n } n instead of 
{A„}„ do not follow from (|3.ip . (|3.4|) . and (|3.5p with {A n } n , but they have to 
be proved independently. See Proposition 3 in [T] for the analogous claim in the 

Case Of f pre , fpost > 0. 

Proposition 3.6 Let {B n } n be an uniformly bounded matrix sequence, with B n 
Hermitian positive definite matrices. For any no > 0, let V^, pre = I n ~ u pre B n , 
Ki,post = In ~ ^postB n be the pre-smoothing and post-smoothing iteration ma- 
trices, respectively considered in the TGM algorithm. Then, there exist as, pre, 
as. post > independent of no such that for all x G C^™ * 1 

||Vo,prc2:||i < |M|| - OB.pjrellVo.prefflUg, (3.12) 

||Vo,posta;||| < \\x\\ 2 Bo - a B ,post\\x\\ B 2. (3.13) 
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Proof. Relation (|3.13p is equivalent to the existence of an absolute positive 
constant as !pos t such that 

(Jo — uj pos tB ) 2 B < B — as, post -Bo, 

i.e., 

■^post^O — 2<x>p OS t-Zo < —OiB,postIo- 

The latter is equivalent to require that the inequality QB lPO st < <^post(2 — w p0 stA) 
is satisfied for any eigenvalue A of the Hermitian matrix Bq with a_B jPOS t > 
independent of uq. Now, let [to, M] be any interval containing the topological 
closure of the union over all n of all the eigenvalues of B n . Then it is enough 
to consider 

as, post < Wpost min (2 - uJ post \) = w post (2 - uj post M), 

\e[m,M] 

where the condition uj pos t < 2/M ensures as !pos t > 0. 

By exploiting an analogous technique, in the case of relation (|3.12[) . it is sufficient 
to consider 

^ . w pro (2 — w pre A) 
a_B,prc < w pre mm — — — 

Ae[o,M] (l-w pre A) 2 



2w prc if < uj pic < 3/(2M), 

if 3/(2M) < Lu prc < 2/M, 



W prc (2 — U! pre M) 



(1 - lo wc MY 

where we consider the only interesting case m = 0, since m > is related to the 
case of well-conditioned systems. • 

In this way, according to the Ruge and Stuben algebraic theory, we have 
proved the TGM optimality, that is its convergence rate independent of the size 
N(n) of the involved algebraic problem. 

Theorem 3.7 Let {B n } n be an uniformly bounded matrix sequence, with B n 
Hermitian positive definite matrices. Under the same assumptions of Proposi- 
tions 13.51 and 13.61 the TGM with only one step of post- smoothing converges to 
the solution of B n x = b and its convergence rate is independent of N(n). 

Proof. By referring to Propositions 13.51 and 13.61 the claim follows according to 
Theorem 13.11 • 
Few remarks are useful in order to understand what happens when also a pre- 
smoothing phase is applied. 

A) The first observation is that the convergence analysis can be reduced some- 
how to the case of only post-smoothing. Indeed, looking at relation (|2.2p 
and recalling that the spectra of AB and BA are the same for any pair 
(A, B) of square matrices (see [7]), it is evident that 

TGM = V^posf CGC Vq° p ^° 

has the same spectrum, and hence the same spectral radius p(-), as 

T/ y 0,pro T /Impost 

K 0,pre K 0,post ^^r^O; 
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where the latter represents a TGM iteration with only post-smoothing. 
Therefore 



p(TGM 



< 1 1 V 0,pi" ^0,posT CGCq 1 1 A 



< 



v/l - «post//3 



where <5 pos t is the post-smoothing constant of the cumulative stationary 
method described by the iteration matrix Vq ° p ^ b Vq ^°^ . 

B) Setting vo, P rc — ^o.post = 1 and with reference to Item A), we easily deduce 

that ap OS t > a pos t where the latter is the post-smoothing constant related 
to the sole post-smoothing method with iteration matrix Vo iPOS t- Further- 
more, if the two iteration matrices Vo jPre and Vo jPOS t are chosen carefully, 
i.e., by taking into account the spectral complementarity principle, then 
we can expect that a pos t is sensibly larger than a p0 st, so that the TGM 
with both pre-smoothing and post-smoothing is sensibly faster than that 
with only post-smoothing. 

C) Items A) and B) show that the TGM iteration with both pre-smoothing 

and post-smoothing is never worse than the TGM iteration with only 
post-smoothing. Therefore Theorem 13. 71 implies that the TGM with both 
post-smoothing and pre-smoothing is optimal for systems with matrices 
B n under the same assumptions as in Theorem 13. 71 

D) At this point the natural question arises: is it possible to handle directly 

assumption (|3.6|) . instead of assumption (|3.2j) ? As observed in Remark l3~4l 
these two assumptions are tied up and indeed they represent the approx- 
imation property on the range of CGCq and unconditional, respectively. 
However, from a technical viewpoint, they are very different and in fact 
we are unable to state a formal analog of Proposition 13.51 by using 
More precisely, for concluding that 



\CGC x\\ Ao <(3 A \\x\\% 



implies 



\CGC x\\ 2 Bo < f3 B \\x\\%. 



with Xo = I as in Remark 13. 31 and with i9, f3 A , Pb absolute constants, and 
A n < ■QBn, we would need X < Y, X, Y > implies X 2 < ^Y 2 with some 
7 positive and independent of n. The latter with 7 = 1 is the operator 
monotonicity of the map z — > z 2 which is known to be false in general 7! . 
We should acknowledge that there exist important subclasses of matrices 
for which X < Y, X, Y > implies X 2 < -fY 2 . However, this matrix 
theoretic analysis of intrinsic interest goes a bit far beyond the scope of 
the present paper and will be the subject of future investigations. 

E) Remark 13.41 furnishes an interesting degree of freedom that could be ex- 
ploited. For instance if we choose Xo — Aq, by assuming suitable order 
relations between {A n }„ and {B n } n , then proving that 

\\CGC x\\ 2 Ao < p A \\x\\ Ao 
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implies 

\\CGC Q x\\% < B \\x\\% 

with Pa, Pb absolute constants, becomes easier, but, conversely, the 
study of the pre-smoothing and post-smoothing properties becomes more 
involved. 

3.2 MGM convergence and optimality: a discussion 

In this section we briefly discuss the same question as before, but in connection 
with the MGM. First of all, we expect that a more severe assumption between 
{A n }„ and {B n } n has to be fulfilled in order to infer the MGM optimality for 
{B n } n starting from the MGM optimality for {A n } n . The reason is that the 
TGM is just a special instance of the MGM when setting 1 = 1. 
In the TGM setting we have assumed a one side ordering relation: here the 
most natural step is to consider a two side ordering relation, that is to as- 
sume that there exist positive constants i?i,i?2 independent of n such that 
$xB n < A n < i?2B n , for every n large enough. The above relationships simply 
represent the spectral equivalence condition for sequences of Hcrmitian posi- 
tive definite matrices. In the context of the preconditioned conjugate gradient 
method (see [5]), it is well known that if {P n }n is a given sequence of optimal 
(i.e., spectrally equivalent) preconditioners for {A„}„, then {P n } n is also a se- 
quence of optimal preconditioners for {B n } n (see e.g. [H]). The latter fact 
just follows from the observation that the spectral equivalence is an equivalence 
relation and hence is transitive. 

In summary, we have enough heuristic motivations in order to conjecture that 
the spectral equivalence is the correct and needed assumption and, in reality, 
the numerical experiments reported in Section [5] give a support to the latter 
statement. 

From a theoretical point of view, as done for the TGM, we start from the Ruge- 
Stiiben tools [25j in the slightly modified version contained in Theorem 2.3 in 
[2] , that is taking into account Remark 13.31 and, for the sake of simplicity, we 
assume no pre-smoothing i.e., v prc = 0. The matrix inequalities coming from 
the assumption (2.9) in [2] are very intricate since they involve simultaneously 
projector operators and smoothers: whence, it is customary to split it into 
the smoothing property (relation (2.11a) in [2]) and the approximation property 
(relation (2.11b) in [5]). As usual the smoothing property does not pose any 
problem. However, we encounter a serious technical difficulty in the second in- 
equality, i.e., when dealing with the approximation property. More precisely, we 
arrive to compare two Hermitian projectors, depending on the same p s s +1 with 
the first involving A ns and the second involving B ria . Unfortunately, they can 
be compared only in very special and too restricted cases: the needed assump- 
tion would not involve ordering, but only the fact that the columns of Al/ s 2 p s s +1 
and those of Bn{ 2 p s s +1 span the same space. 

However, as already mentioned, the numerical tests tell us that the latter diffi- 
culty is only a technicality and that the right assumption should involve spectral 
equivalence. Therefore, in future investigations, other directions and proof tech- 
niques have to be explored. 
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4 A case study: discrete Laplacian-plus-diagonal 

systems 

In the present section, we analyze a specific application of the results in Section 
[31 More precisely, we consider a multigrid strategy for solving Laplacian-plus- 
diagonal linear systems arising from standard Finite Differences (FD) discretiza- 
tions of the problem 

Lu = -V 2 m(x) + fj,(x)u(x) = h(x), x e Q,, (4.1) 

where n{x) and h{x) are given bounded functions, f2 = (0, l) d , d > 1, and with 
Dirichlet, periodic or reflective boundary conditions. Thus, we are facing with 
a matrix sequence 

{B n } n = {A n + D n } ni (4.2) 

where the structure of the matrix sequence {j4„}„ is related both to the FD 
discretization and to the type of the boundary conditions and where {D n } n is 
a sequence of uniformly bounded diagonal matrices, due to the hypothesis that 
n(x) is bounded. 

Since a fast TGM and MGM working for the Toeplitz (or r - the r class is 
the algebra associated to the most known sine transform [8]) part is well-known 
[T51 [TBI [U [57] , we are in position to apply the tools in the preceding section 
in order to show that the same technique works, and with a cost linear as the 
dimension, in the context (|4.2p too. In the same way, the extension of suitable 
MGM procedures proposed in the case of the circulant [30], DCT-III cosine 
[T21 I3"5] , or t [ini HH] algebra, can be considered according to the corresponding 
boundary conditions. Clearly, this case study is just an example relevant in 
applications, while the results in Section [3] are of much wider generality. 
Once more, we want to remark that, unfortunately, there is a gap in the theory 
with regard to the MGM, even if the numerical tests reported in Section [5] 
suggest that the MGM applied to matrices in {B n — A n + D n } n is optimal 
under the assumptions that the same MGM is optimal for {A n } n , A n symmetric 
positive definite matrix, and {D n } n uniformly bounded matrix sequence, with 
A n < "&B n uniformly with respect to n and with some fixed > independent 
of n. Clearly if the matrices D n are also nonnegative definite then the constant 
■d can be set to 1. This result can be plainly extended to the case in which D n 
is a (multilevel) banded correction. 

4.1 One-Dimensional case 

According to the FD approximation of (|4.1|) with Dirichlet boundary conditions, 
we obtain the matrix sequence 

{B n } n = {A n + D n } n 

where {A„}„ = {tridiag„ [— 1,2, — l]} n and {D n }„ is a sequence of diagonal 

(n) 

matrices whose diagonal entries d\ , i — l,...,n, are uniformly bounded in 
modulus by a constant M independent of n. Since 



11 



we impose the condition 

n mm a) + n > c 

l<i<n 1 

for some c > independent of n (we consider only the case mini<j< n < 0, 
since the other is trivial). Thus, also {B n } n is an uniformly bounded positive 
definite matrix sequence and 

2 

so satisfying the crucial assumption A n < "t}B n in Proposition l3.5l with d = ir 2 /c. 
Let us consider B € t' loX " (l , with 1-index n > 0. Following [13 [27], we denote 
by To 1 e E" oXni , n Q = 2n x + 1, the operator such that 

(T u _ / 1 for i = 2 3, i = l,...,ni, , > 

I^oKj - | o otherwise, 1 ' 

and we define a projector (p\) H , p\ € R n ° Xni as 

pj = ^oT \ Po = tridiag [1, 2, 1]. (4.4) 

Thus, the basic step in order to prove the TGM optimality result is reported in 
the proposition below. It is worth stressing that the claim refers to a tridiagonal 
matrix correction, since, under the quoted assumption, each diagonal correction 
is projected at the first coarse level into a tridiagonal correction, while the 
tridiagonal structure is kept unaltered in all the subsequent levels. 

Proposition 4.1 Let Bq = tridiag [-1,2,-1]+T £ R n ° XIl ° ; with n a > and 
Tq being a symmetric uniformly bounded tridiagonal matrix such that Aq < dBo, 
with Aq = tridiag [—1,2,-1] and some d > 0. Letpl = (l/-v/2)tridiag [1,2, 1]T X , 
with no = 2n\ + 1. Then, 

(pl) H B Po = tridia gl [-1,2, -1] + T t 

where T\ £ l" 15 " 11 is a symmetric uniformly bounded tridiagonal matrix with 
Ax < #Bx, Pi = {pl) H B pl, A 1 = tridia gl [-1,2, -1]. 

Proof. For the Toeplitz part refer to [57] ■ For the tridiagonal part we just 
need a simple check. In fact, the product PoPoPo is a 7-diagonal matrix (Po 
and To are tridiagonal) and the action of Tq, on the left and on the right, 
selects the even rows and columns so that the resulting matrix is still tridi- 
agonal. Since Aq < i}B , Aq = tridiag [—1,2,-1], B\ = (pJ) H P Po, and 
Ai = tridiag-L [-1,2,-1] = (pl) H AqpI it is evident that A x < $B 1 . Finally, 
the uniform boundedness is guaranteed by the uniform boundedness of all the 
involved matrices. • 

Corollary 4.2 Let B Q = tridiag [-1,2,-1]+T € W l ° xn ° , withn Q > andT 
symmetric tridiagonal matrix such that Aq < "&Bq, with Aq = tridiag [—1,2,-1] 
and some d > 0. Let pi = (l/^tridiagQ [1, 2, IJTq 1 , with n Q = 2n x + 1. Then, 
there exists /3b > independent of uq so that inequality (13. 2|) holds true. 
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Proof. Let A = tridiag [-1,2,-1] € M"» x "». Then relation is fulfilled 
with the operator defined in (|4.4[) and with a certain /J^ independent of no, 
as proved in [27]. Moreover, from the assumption we have Ao < {)Bq so that 
Proposition 13.51 implies that (|3.2|) holds true for B with a constant (3b = fipA- 



Corollary 4.3 Let {B n } n be the sequence such that B n = tridiag n [—1, 2, — 1] + 
T„ G R nxn with T n symmetric uniformly bounded tridiagonal matrices. Then, 
there exist as iPrc , as jPOS t > independent of n, so that inequalities (|3.12jl and 
(jSIIll) hold true. 

Proof. It is evident that {B n } n is a sequence of symmetric positive definite 
matrices uniformly bounded by 4 + M, with ||T n ||2 < M independent of n, so 
that the thesis follows by the direct application of Proposition 13. 61 • 



4.2 Two-Dimensional case 

Hereafter, we want to consider the TGM and MGM extension to the case 
d > 1. Due to the discretization process, it is natural, and easier, to work with 
d— indices n — (n^\ . . . , nW), with integer positive number, r = 1, . . . ,d. 

In this case the matrix dimension is iV(no) = Yit=i n o an<4 wnen considering 
the projected matrices of size N(ni) we have that ni is again a d— index and 
we assume not only N(nx) < N(n ), but also m < n componentwise. 
We discuss in detail the two-level case, since the d— level one is a simple gener- 
alization. Thus, in the two-level case, we are dealing with the matrix sequence 

{B n } n = {A n + D rL \ n 

where {A n } n = {tridiag^i) [-1,2, -l]®/ n( 2) +I n a) ® tridiag„ ( 2) [-1,2, -l]} n 
and {D n } n is a sequence of diagonal matrices whose diagonal entries d\ n \ i — 
1, . . . , N(n), are uniformly bounded in modulus by a constant M independent 
of n. Since 

X 2 



A m in(A) = 4 sin rT-TTT-TT + 4 sin 



n 



2(n( 1 ) + l)y' \2(n( 2 ) + l) 



W]2 [n(!)] 2 



we impose the condition 



ib 2 min d-™' + 7r 2 > c 

l<i<N(n) 



for some c > independent of n, with ^ = mia^{7i^}. Thus, also {B n } n is an 
uniformly bounded positive definite matrix sequence and 



c — 7T 2 c 



so satisfying the crucial assumption A n < dB n in Proposition l3.5l with d = ir 2 /c. 
The projector definition can be handled in a natural manner by using tensorial 



13 



arguments: {pq) h is constructed in such a way that 
Po = Port 

P Q = tridiag m [1,2,1]® tridiag (2) [1,2,1], 

rt = TZ(nV)®Tt(n™) 

with — 2n{^ + 1 and where Tq 1 (n^ ) £ R"o ' x "i ' is the unilevel matrix given 
in (|4.3p . Notice that this is the most trivial extension of the unilevel projector 
to the two-level setting and such a choice is also the less expensive from a 
computational point of view: in fact, p^ = tq((2 + 2cos(£i)(2 + 2cos(i2)))^d 
equals [r^, (p(2 + 2 cos(t 1 )))T 1 (4 1) )] ® [r n ( 2 , (p(2 + 2 cos^)))^ 1 ^)]. 
The proposition below refers to a two-level tridiagonal correction for the same 
reasons as the unilevel case. 

Proposition 4.4 Let 

B Q = tridiagm [-1,2,-1]®/ (2 >+7 (1) ® tridiagm [-1,2,-1]+T £ K ^(™o)xJV(„ o); 

with no > and To being a symmetric uniformly bounded tridiagonal block ma- 
trix with tridiagonal blocks such that Ao < $Bo, with Ao = tridiag„(i) [—1, 2, —1]® 
I n (2) + ® tridiag„(2) [—1,2,-1] and some ■& > 0. Let 

p l = (tridiag^, [1,2,1]® tridiag n ( 2) [1, 2, l]){T^n^) ® ijftng )), 

with n Q r ^ = 2n r ^ + 1, r = 1, 2. Then, B\ :— (pq) h Bop\ coincides with A\ + T\ 
where T\ £ R N ( n i) xN ( n i) j s a symmetric uniformly bounded tridiagonal block 
matrix with tridiagonal blocks and where A\ is a two-level r tridiagonal block 
matrix with tridiagonal blocks asymptotic to 

tridiag n (i) [-1, 2, -1] ® 1^) + J^i) ® tridiag n m [-1,2,-1], 
so i/iai Ai < . 

Proof. For the r part refer to [27] • For the two- level banded part it is a simple 
check. In fact, the product PoToPo is a 7-diagonal block matrix with 7-diagonal 
blocks (Po is a tridiagonal block matrix with tridiagonal blocks) and the action 
of J7q, on the left and on the right, selects even rows and columns in even 
blocks with respect to the rows and columns, so that the resulting matrix has a 
tridiagonal block pattern with tridiagonal blocks. The order relation follows as 
a direct consequence of the assumption Aq < $Bo and the uniform boundedness 
is implied by the uniform boundedness of all the involved matrices. • 

Corollary 4.5 Let B =A +T £ M^(™o)xjv(n ) with nQ > 0) 

Ao = tridiag m [-1, 2, -1] ® / w + 1 m ® tridiag (2) [-1,2,-1], 

and Tq symmetric tridiagonal block matrix with tridiagonal blocks such that Aq < 
•QBo for some § > 0. Let 

pi = (tridiag n (i, [1,2,1]® tridiag n ( 2) [1, 2, l]){T^n^) ® 

(r) (r) 

with n = 2n\ + 1, r = 1,2. Then, there exists /3b > independent of no so 
that inequality (|3.2p holds true. 
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Proof. The proof can be done following the same steps as in Corollary 14. 21 • 

Corollary 4.6 Let {B„} n be the sequence such that B n = A n + T n with 

A n = tridiag n( i) [-1,2,-1] ® I nW + I n m ® tridiag„( 2) [-1,2,-1], 

and with T n symmetric uniformly bounded tridiagonal block matrix with tridi- 
agonal blocks. Then, there exist a_B jPrc , aB,p re > independent of n, so that 
inequalities (I3.12p and (|3.13[) hold true. 

Proof. The proof can be worked out as in Corollary 14.31 since the sequence 
{B n } n is uniformly bounded by 8 + M, with ||T„|| 2 < M independent of n by 
assumption. • 



5 Numerical Examples 

We test our TGM and MGM (standard V-cycle according to Section [5]) for 
several examples of matrix corrections {£>„}„, D n £ ^N{n)xN(n) ^ N(n) = 

We will consider nonnegative definite band corrections and indefinite band 
corrections. By referring to Section HJ the case of nonnegative definite correc- 
tions implies trivially that A n < B n so that the desired constant is i9 = 1. 
However, as observed in real- world applications (see [14]). the most challenging 
situation is the one of indefinite corrections. 

Concerning nonnegative definite corrections, the reference set is defined accord- 
ing to the following notation, in the unilevel and in the two-level setting, re- 
spectively: 











dO 


dl 


d2 


d3 


d4 


ID 


s = 


1,.. 


.,N{n) 





s 

s+1 


|sin(s)| 


|sm(s)|f 3T i 


S 

N{n) 




s = 


(<- 


i)m + j 





i 

i+1 


|sin(i)| 


IsiuWI^i 


s 

N{n) 


2D 


i = 


1,.. 


. ,ni 




+ I+T 


+ |sin(j)| 


+ |sin(j)|^i 






i = 


1,.. 


-,n 2 













The case of indefinite corrections is considered in connection with Laplacian 
systems with Dirichlet boundary conditions: in that setting the diagonal entries 
di of D n are generated randomly. Finally higher order differential operators 
and linear systems arising from integral equations in image restoration are con- 
sidered at the end of the section. 

The aim is to give numerical evidences of the theoretical optimality results of 
TGM convergence and also to their extension in the case of the MGM applica- 
tion. 
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The projectors are properly chosen according to the nature of structured part, 
while we will use, in general, the Richardson smoothing/intermediate iteration 
step twice in each iteration, before and after the coarse grid correction, with 
different values of the parameter lo. 

According to the definition, when considering the TGM, the exact solution of 
the system is obtained by using a direct solver in the immediately subsequent 
coarse grid dimension, while, when considering the MGM, the exact solution of 
the system is computed by the same direct solver, when the coarse grid dimen- 
sion equals 16 d (where d = 1 for the unilevel case and d — 2 for the two-level 
case) . 

In all tables we report the numbers of iterations required for the TGM or MGM 
convergence, assumed to be reached when the Euclidean norm of the relative 
residual becomes less than 10~ 7 . We point out that the CPU times are consis- 
tent with the iteration counts. 

Finally, we stress that the matrices A n at every level (except for the coarsest) 
are never formed since we need only to store the nonzero Fourier coefficients of 
the generating function at every level for matrix-vector multiplications. Thus, 
besides the 0(N(n)) operations complexity of the proposed MGM both with 
respect to the structured part and clearly with respect to the non-structured 
one, the memory requirements of the structured part are also very low since 
there are only O(l) nonzero Fourier coefficients of the generating function at 
every level. On the other hand, the projections of the initial diagonal correction 
are stored at each level according to standard sparse matrix techniques during 
the pre-computing phase. 

5.1 Discrete Laplacian-plus-diagonal systems 

The numerical tests below refer to convergence results in the case of matrix 
sequences arising from the Laplacian discretization, in the unilevel and in the 
two-level settings, respectively. 

5.1.1 Dirichlet boundary conditions 

Firstly, we consider the case of Dirichlet boundary conditions so that the ob- 
tained matrix sequence is the Toeplitz/r matrix sequence {T„(/)} n generated by 
the function f(t) — 2 — 2 cos(i), t £ (0, 2ir]. The projector is defined as in (|4.3|) 
and (|4.4|) . while the parameters lo for the smoothing/intermediate iterations are 
chosen as 

1 , 1 

WprC " 2(11/1100 + HA.IU) Wp ° St " II/IIoc + IIAJoc' 

with f prc = ^ post = 1. 

The results in Table [T] confirm the optimality of the proposed TGM in the sense 
that the number of iterations is uniformly bounded by a constant not depending 
on the size N (n) indicated in the first column. Moreover, it seems that this claim 
can be extended to the MGM convergence. Notice, also, that the number of 
iterations is frequently the best possible since it equals the number of TGM 
iterations. 

The case of the diagonal correction dA deserves special attention: as shown 
in the first column, just one pre-smoothing/intermediate and post-smoothing 



16 



B n = tridiag n [—1,2,-1]+ Diagonal 



TGM 




MGM 


N(n) 
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dl 
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AT(n) 


dO 


dl 


d2 


d3 


d4 


























p = 


p=l 


31 


2 


7 


7 


7 


7 
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7 
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8 


16 
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Table 1: Number of iterations required by TGM and MGM - unilevel cases 
(refer to (|5.ip for the definition of the constant p). 
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10 
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10 
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511 2 


16 


10 


13 


13 
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16 
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12 


12 


36 


16 



Table 2: Number of iterations required by TGM and MGM - two- level cases 
(refer to (|5.1j) for the definition of the constant p). 



iteration at each coarse grid level are not sufficient to ensure the optimality. 
Moreover, it is enough to consider a trick that keeps unaltered the 0(N(n)) 
computational cost as proved in [30j (only the multiplicative constant hidden 
in the big O can increase): at each projection on a coarser grid the number of 
smoothing iterations performed at that level is increased by a fixed constant p, 
i.e., according to the MGM notation of Section [21 we set 

v s+ i = v s + p, s = 0, . . . , I - 1, vq = 1. (5.1) 

The optimality result in the second column relative to MGM in the c?4 case 
is obtained just by considering p = 1. This phenomenon is probably due to 
some inefficiency in considering the approximation ||D n ||oo m the tuning of the 
parameter w pro and w pos t- In fact, it is enough to substitute, for instance, the 
post-smoother with the Gauss-Seidel method in order to preserve the optimality 
also for p = 0. 

Other examples of Toeplitz/r linear systems plus diagonal correction can be 
found in |23j , corresponding to Sinc-Galerkin discretization of differential prob- 
lems according to [20] • 

By using tensor arguments, our results plainly extend to the two-level setting 
and the comments concerning Table [2] are substantially equivalent as in the 
unilevel case. 

Before dealing with other type of boundary conditions, we want to give a 
comparison of the performances of the proposed method with respect to those 
achieved by considering, for instance, the conjugate gradient (eg) method. Table 
131 reports, for increasing dimension, the Euclidean matrix condition number 
&2(^4n + D n ), together with the number of iterations required by the eg. As 
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B„ = tridiag n 
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Table 3: Euclidean condition number ki{A n + D n ) and number of iterations 
required by eg - unilevel and two-level cases. 



well known in the case dO, the eg method requires all the N(n) steps in order 
to reach the convergence. Moreover, the non-structured part in the cases dl, 
d2, d3 increases the minimum eigenvalue of the resulting matrix so that the 
whole condition number becomes moderate. As a consequence the standard eg 
method is also effective. Notice that this good behavior is no longer observed 
in the case d4, while our MGM technique is still optimal. The same trend is 
observed in the two-level setting. 



5.1.2 Randomly generated indefinite corrections 

As a further interesting case, we want to test our proposal in the case of ran- 
domly generated matrix corrections. More specifically, we consider diagonal, 
symmetric tridiagonal, symmetric pentadiagonal matrix corrections with ran- 
dom entries uniformly distributed on the unit interval (cases d5, d7, and d9, 
respectively) or normally distributed with mean zero and standard deviation 
one (cases d6, d8, and dlO, respectively). Notice that in such a way we are also 
considering indefinite corrections. Thus, in order to obtain a positive definite 
matrix B n and in order to satisfy the crucial relation A n < d-B„ for some pos- 
itive d independent of n, the arising random matrices corrections are suitable 
scaled by l/^n 2 ) in the unilevel setting, with 7 being the number of non-zero 
diagonals, and by l/^nA 1 )) 2 ) in the two-level setting (assuming = n^). 
Table 0] reports the Euclidean condition number and the mean of the number 
of iterations required by the MGM in the unilevel and two-level setting by con- 
sidering, for each case, ten examples of random matrix corrections. 
All these results confirm the effectiveness of our proposal. Though the Eu- 
clidean condition numbers are fully comparable with those of the dO case, the 
number of required iterations does not worsen. Conversely, the eg method re- 
quires for instance in the d5 case N(n) iterations in the unilevel setting, and 
83, 163, 318, 621, 1212 in the two-level one. 

It is worth stressing that the pentadiagonal corrections are reduced at the first 
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B n = tridiag n [— 1, 2, — l]+random correction 
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B n = tridiag n( i) [-1, 2, -1] <gl I n{2) + I n (i) ® tridiag n(2 ) [— 1, 2, -l]+random correction 



N(ri) 


d5 




d6 






d7 




d8 




d9 




dlO 






k 2 


nit 


k 2 


nit 


k 2 




nit 


k 2 


nit 


k 2 


nit 


k 2 


nit 


31 2 


4.02e+2 


16 


4.15c+2 


16 


3.98e 


f2 


16 


4.150+2 


16 


3.96e+2 


16 


4.14e+2 


16 


63 2 


1.61e+3 


16 


1.66e+3 


16 


1.59c 


f3 


16 


1.66e+3 


16 


1.59e+3 


16 


1.65e+3 


16 


127 2 


6.47e+3 


16 


6.64e+3 


16 


6.39e 


f3 


16 


6.63e+3 


16 


6.36e+3 


16 


6.63e+3 


16 


255 2 


2.58e+4 


16 


2.65e+4 


16 


2.55e 


+-4 


16 


2.65e+4 


16 


2.54e+4 


16 


2.65e+4 


16 


511 2 


1.03e+5 


16 


1.06e+5 


16 


1.02c 


+5 


16 


1.06e+5 


16 


1.01e+5 


16 


1.06e+5 


16 



Table 4: Euclidean condition number k2(A n + D n ) and mean number of itera- 
tions required by MGM - unilevel and two-level cases. 



projection to tridiagonal matrices. More in general, bigger patterns are reduced 
after few steps to a fixed pattern driven by the projector pattern (see [UJ [TJ [3] ) . 



5.1.3 Periodic and Reflective boundary conditions 

In the case of periodic boundary conditions the obtained matrix sequence is 
the circulant matrix sequence {S n (f)} n generated by the function f(t) = 2 — 
2cos(t),i G (0,27t]. Following [30], we consider the operator Tq 1 e R n » xni , 
tiq = 2rt\, such that 

(T i^ f 1 for i = 2j - 1, j = 1, . . . ,m, 
1 oKj \ otherwise, 

and we define a projector (p\) H , € M n ° xni , as 

pi = PoT* , P = So (?) , P (t) = 2 + 2 cos(t) . 

It must be outlined that in the dO case the arising matrices are singular, so that 
we consider the classical Strang correction [33] 

Sn (f) = s no (f) + / (jv^y) 

where e is the vector of all ones. The results in the top part of Table [5] confirm 
the optimality of the proposed TGM and its extension to MGM (the case dA 
requires to set p = 4). 

When dealing with reflective boundary conditions, the obtained matrix sequence 
is the DCT III matrix sequence C n (f) n generated by the function f(t) = 2 — 
2cos(t),i € (0,27r]. Following [12], we consider the operator T X € r° xni , 
no = 2ni, such that 

(T u _ / 1 for ie{2j -l,2j}, j = l,...,m, 
[ oJjJ \ otherwise, 
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B n = unilevel circulant S n (f)+ Diagonal, f(t) = 2 — 2cos(i) 



TGM 




MGM 


N(n) 


dO 


dl 


d2 


d3 


d4 




AT(n) 


dO 


dl 


d2 


d3 


d4 


























p = 


p = 4 


32 


2 


6 


7 


7 


7 




32 


2 


6 


7 


7 


7 


7 


64 


2 


6 


7 


7 


7 




64 


7 


6 


7 


7 


10 


7 


128 


2 


6 


7 


7 


7 




128 


7 


6 


7 


7 


16 


7 


256 


2 


6 


7 


7 


7 




256 


8 


6 


7 


7 


22 


7 


512 


2 


6 


6 


7 


6 




512 


8 


6 


6 


7 


29 


8 






B n = 


unilevel DCT III C„(/)+ Diagonal, /(t) = 


= 2 - 2cos(t) 






























TGM 




MGM 


N(n) 


dO 


dl 


d2 


d3 


d4 




N(n) 


dO 


dl 


d2 


d3 


d4 


























p = 


p = 2 


32 


7 


6 


7 


7 


6 




32 


7 


6 


7 


7 


6 


6 


64 


7 


5 


6 


6 


5 




64 


7 


5 


6 


6 


7 


5 


128 


7 


5 


7 


7 


5 




128 


7 


5 


7 


7 


11 


5 


256 


7 


4 


7 


7 


4 




256 


7 


5 


6 


6 


17 


5 


512 


7 


4 


6 


6 


4 




512 


7 


4 


7 


7 


27 


6 



Table 5: Number of iterations required by TGM and MGM - unilevel cases 
(refer to (|5.ip for the definition of the constant p). 



and we define a projector (pg) H , Pq G 

PoT \ P = C (p), P (t) 



as 



Pi 



2cos(£). 



Again, the results in bottom part of Table confirm the optimality of the 
proposed TGM and its extension to MGM. It is worth stressing that in the dO 
case we are considering the matrix 



Cn (f) ~ C m {f) + / ) 



Furthermore, the case d4 requires to set p = 2 in order to observe optimality. 
By using tensor arguments, our results plainly extend to the twodevel setting 
and the comments concerning Table [6] are substantially equivalent as in the 
corresponding unilevel case. 



5.2 Other examples 

In this section we give numerical evidences of the optimality of TGM and MGM 
results in a more general setting. 



5.2.1 Higher order r discretizations plus diagonal systems 

We consider r matrix sequences arising from the discretization of higher order 
differential problems with proper homogeneous boundary conditions on dfl: 

d fp-q 

i.e., {B n = A n + D n } n , where A n = r n (/) with f(t) = £ti( 2 - 2cos(^)) 9 - 
More specifically, in the unilevel case we define p{t) = [2 + 2 cos(t)\ w where w is 
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B n = two-level circulant S n (f)+ Diagonal, /(ii,i2) = 4 — 2cos(ti) — 2cos(<2) 



TGM 




MGM 


N(n) 


dO 


dl 


d2 


d3 


d4 




N(n) 


dO 


dl 


d2 


d3 


d4 
p = 


p=l 


32 2 


15 


8 


11 


11 


14 






15 


8 


11 


11 


14 


14 


64 2 


15 


7 


11 


11 


15 




64 2 


15 


7 


11 


11 


15 


15 


128 2 


15 


7 


11 


11 


15 




128 2 


15 


7 


11 


11 


16 


14 


256 2 


15 


7 


11 


11 


15 




256 2 


15 


7 


11 


11 


24 


14 


512 2 


15 


7 


11 


11 


15 




512 2 


15 


7 


11 


11 


34 


14 



B n = two-level DCT III C„(/)+ Diagonal, /(ti,t 2 ) = 4 - 2cos(*i) - 2cos(t 2 ) 



TGM 




MGM 


N(n) 


dO 


dl 


d2 


d3 


d4 




N(n) 


dO 


dl 


d2 


d3 


d4 


























p = 


p=l 


32 2 


16 


6 


10 


10 


12 




32 2 


16 


6 


10 


10 


12 


12 


64 2 


16 


6 


10 


10 


11 




64 2 


16 


6 


10 


10 


11 


11 


128 2 


16 


5 


10 


10 


11 




128 2 


16 


5 


10 


10 


11 


10 


256 2 


16 


5 


9 


9 


11 




256 2 


16 


5 


9 


9 


17 


9 


512 2 


16 


5 


9 


9 


11 




512 2 


16 


5 


9 


9 


27 


9 



Table 6: Number of iterations required by TGM and MGM - two-level cases 
(refer to (|5.ip for the definition of the constant p). 



chosen according to conditions in [T51 [HI [27] : in order to have a MGM optimality 
we must take w at least equal to 1 if q = 1 and at least equal to 2 if q = 2, 3. 
Clearly, the lower is the value of w, the greater will be the advantage from a 
computational viewpoint. Indeed, Table [7] confirms the need of these constraints 
with respect to the case q — 2, this being the only dO case where we observe 
a growth in the number of iterations with respect to N(n). Nevertheless, it 
should be noticed that in the same case the contribution of the non-structured 
part improves the numerical behavior since the minimal eigenvalue is increased. 
The remaining results in Table [7J confirm the optimality of the corresponding 
MGM (the dA case requires to set p in a proper way as just observed in the 
Laplacian case). 

Notice that the bandwidth of the non-structured diagonal correction is increased 
by subsequent projections until a maximal value corresponding to Aw — 1 is 
reached (for a discussion on the evolution of the bandwidth when a generic 
(multilevel) band system is encountered see [T31 IT] 12"]). 

With respect to the two-level problem, we consider again the most trivial ex- 
tension (and less expensive from a computational point of view) of the unilcvcl 
projector to the two-level setting, given by P n = r n (p) with p(ti,t 2 ) = [(2 + 
2cos(£i))(2 + 2 cos^))]™, w = 1,2,3. The comments concerning the two-level 
setting in Table [7J are of the same type as in the unilcvcl one. 

5.2.2 Higher order circulant discretizations plus diagonal systems 

We consider circulant matrix sequences arising from the approximation of higher 
order differential problems with proper homogeneous/periodic boundary condi- 
tions on 9f2 as in (|5.2p . i.e., {B n = A n + D n } n , where A n — S n (f) with 
f(t) = 53i=i(2 — 2cos(ii)) 9 . The choice of the generating function for the pro- 
jector is the same as in the previous section (see [30] )• Indeed, Table [5] shows 
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B n =unilevel r+ Diagonal, f(t) = (2 - 2cos(t))« 















q = 


o 
Z 
















w = 1 








w — 


z 








N(n) 


dO 


dl 


d2 


d3 


d4 


AT/*-, \ 

IV [n ) 


dO 


A 1 


An 
dz 


do 




d4 
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p = 4 












p — u 
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31 


20 


13 


16 


16 


13 


13 


31 


16 


13 




14 


13 




13 


63 


45 


13 


17 


17 


17 


13 


63 


16 


12 


14 


14 


12 




12 


127 


84 


12 


15 


15 


26 


13 


127 


16 


12 


13 


13 


13 




12 


255 


149 


12 


16 


16 


36 


13 


255 


16 


12 


13 


13 


27 




12 


511 


253 


12 


16 


16 


42 


13 


511 


16 


11 


13 


13 


35 




13 














q = 


3 
















w = 2 








w = 


3 








N(n) 


dO 


dl 


d2 


d3 


d4 


JV(n) 


dO 


dl 


d2 


d3 




d4 












p = 


p=l 












p=0 




p=l 


31 


35 


31 


33 


33 


32 


32 


31 


34 


31 


33 


33 


32 




32 


63 


35 


31 


33 


33 


31 


31 


63 


34 


31 


33 


33 


31 




31 


127 


35 


31 


32 


32 


31 


32 


127 


34 


30 


32 


32 


31 




31 


255 


35 


30 


32 


32 


37 


31 


255 


34 


30 


32 


32 


37 




30 


511 


35 


30 


32 


32 


39 


31 


511 


34 


29 


31 


32 


42 




30 



B n =two-level r+ Diagonal, f(t lt t 2 ) = (2 - 2cos(ti)) 9 + (2 - 2cos(t 2 )) 9 
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37 


29 


32 


32 


36 


36 


31 2 


35 


28 


31 


31 


34 


34 


63 2 


44 


28 


32 


32 


43 


35 


63 2 


36 


28 


31 


31 


34 


34 


127 2 


80 


28 


31 


31 


73 


35 


127 2 


36 


27 


30 


30 


56 


34 


255 2 


140 


27 


30 


30 


109 


35 


255 2 


36 


27 


30 


30 


89 


33 


511 2 


235 
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30 


30 


151 


35 


511 2 


36 


27 


30 


30 


129 


33 
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72 


64 


67 


67 


69 


69 


31 2 


68 


61 


64 


64 


66 


66 


63 2 


73 


65 


68 


68 


71 


71 


63 2 


72 


64 


67 


67 


70 


70 


127 2 


73 


64 


68 


68 


137 


71 


127 2 


72 


63 


67 


67 


127 


70 


255 2 


73 


64 


67 


67 


193 


70 


255 2 


72 


63 


67 


67 


182 


70 


511 2 


73 


64 


67 


67 


276 


70 


511 2 


72 


63 


67 


67 


260 


70 



Table 7: Number of required MGM iterations - unilevel and two-level cases 
(refer to (|5.1j) for the definition of the constant p). 



the importance of these constraints with respect to the case dO with q = 2. It 
is worth mentioning that the optimality of the corresponding MGM is again 
confirmed (for the case d4 the parameter p has to be set in a proper way). The 
comments concerning the two- level setting in Table [5] are of the same type as in 
the unilevel one. 

5.2.3 Reflective BCs discretizations plus diagonal systems 

We consider an example of DCT-III matrix sequences arising from the dis- 
cretization of integral problems with reflective boundary conditions (see [22j). 
i.e., {B n — A n + D n } n , where A n = C n (f) with / having nonnegative Fourier 
coefficients as it is required for the point spread function in the modeling of 
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B n = unilevel circulant + Diagonal, fit) = (2 — 2cos(t)) 9 















<? = 


o 
Z 














w = 1 
















N(n) 


dO 


dl 


d2 


d3 


d4 


AT/*-, \ 

IV [n ) 


dO 


A 1 




A Q 

d3 


d4 












p = 


p = 4 












p = 


p = 4 


32 


19 


12 


13 


14 


15 


15 


32 


15 


12 


13 


13 


14 


14 


64 


41 


12 


14 


15 


29 


17 


64 


15 


12 


13 


13 


17 


14 


128 


77 


12 


14 


14 


47 


19 


128 


15 


12 


13 


13 


34 


14 


256 


137 


12 


14 


14 
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29 


30 


30 
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31 
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30 
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30 


30 
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31 
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28 
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B n = two-level circulant + Diagonal, f{t\,ti) = (2 - 2cos(ti)) 9 + (2 - 2cos(t 2 )) 9 
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25 


26 
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33 
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33 


22 


25 
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31 
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22 


25 


25 
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33 


22 


25 


25 
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67 


54 
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66 


53 


57 


58 
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65 


256 2 
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63 
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67 
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57 


57 
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66 
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56 


57 


149 
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Table 8: Number of required MGM iterations - unilevel and two-level cases 
(refer to (|5.1j) for the definition of the constant p). 



image blurring, see [B]. A simple model is represented by f(t) = fd(t) := 
X)i=i(2 + 2cos(fj)) where, by the way, the product /i(£i)/i(t2) is encountered 
when treating super- resolution or high resolution problems, see e.g. |21j . The 
choice of the generating function for the projector is the same as in [12] . 
The results in Table [5] confirm again the optimality of the corresponding MGM 
(the case e?4 requires to set p in a proper way). The observations regarding the 
two-level setting are in the same spirit as those of the unilevel one. 
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B n = unilevel DCT III C„(/)+ Diagonal, f(t) = 2 + 2cos(t). 
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10 


4 
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9 


30 
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two-level DCT III C„(/)+ Diagonal, /(ti,t 2 ) = (2 + 2 cos(*i )) + (2 - 2cos(t 2 )). 



TGM 
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9 


6 


128 2 


4 
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4 
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9 


9 
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5 
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9 
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Table 9: Number of iterations required by TGM and MGM - unilevel and two 
level cases (refer to (|5.ip for the definition of the constant p) . 



6 Concluding Remarks 

The algebraic tools given in Section [3] and Section |4] revealed that, if a suitable 
TGM for a Hermitian positive definite matrix sequence {^4„}„ is available and 
another Hermitian positive definite uniformly bounded sequence {B n } n is given 
such that A n < $B n , for n large enough, then the same strategy works almost 
unchanged for {B n } n too. As an example, this means that if the method is 
optimal for the first sequence then it is optimal for the second as well. The 
same results should hold for the MGM procedures, but here only a wide set 
of numerical evidences has been provided for supporting the claim: the related 
theory will be a subject of future investigations taking into account the final 
remarks in Section 13.11 and the discussion in Section 13.21 

We point out that the latter goal is quite important. Indeed, it is not difficult to 
prove relations of the form d\A n < B n < d2A n with B n being discretization of 
an elliptic variable coefficient problem, A n being the same discretization in the 
constant coefficient case, and where $i,t?2 are positive constants independent 
of n and mainly depending on the ellipticity parameters of the problem. There- 
fore, the above mentioned results would represent a link for inferring MGM 
optimality on a general (possibly high order) variable coefficient elliptic prob- 
lem, starting from the MGM optimality for the structured part, i.e., the one 
related to the constant coefficient discretization. 

Finally, we point out that the latter idea has been used essentially for the struc- 
tured plus diagonal systems coming from approximated elliptic partial differen- 
tial equations with different boundary conditions. However, the same approach 
is applicable to a wide variety of cases, as sketched for instance in Section l5.2.3l 
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A Appendix 



For reader convenience, we report the essential steps of the proof of Theorems 
O and OH] 

Let us start by proving Theorem 13. II As demonstrated in Theorem 5.2 in [25j . 
the existence of (3 > such that 

min \\x-dy\\ 2 Do <l3\\x\\%, Vx e C N ^ 

implies the validity of the so called approximation property only in the range of 
CGCq, i.e., the existence of (3 > such that 

\\CGC x\\ 2 Ao < (3 \\CGC x\\ 2 AoD . lAo Wx e C N ^ (A.l) 
where CGC = I - p\A^ [pl) H A . 

Thus, by virtue of the post- smoothing property (|3.ip and of (|A.1[) . for all x £ 
C^Oo) we find 



\V , post CGC x\\ 2 Ao < \\CGC x\\ Ao ~a post \\CGC x 

< \\CGC xf Ao -^\\CGC xf A 



A D L A 



< [1-^)11*, (A.2) 



1 _ r^i j \\CGC x\\ Ao 



being \\CGC \\ Ao = 1. 

Since TGMq = Vo lPO stGGCo in the case where no pre-smoothing is considered, 



the latter is the same as ||TGAfo|U < y 1 — ap ^ st and hence Theorem 13.11 is 
proved. 

Now let us prove Theorem 13.21 Since the approximation property ()3.6p im- 
plies clearly (| A. 1[) . by repeating the very same steps as before and exploiting 
the post-smoothing property (|3.5p . for all x £ (£ N ( n o) W e find 

\\V , post CGC V , we x\\ Ao < \\CGC Vo, P r C x\\ Ao - a post ||C , GC V& J p ro s||^ oD -i > i 

[ "'" * "CGC V , pre xf Ao . (A.3) 



3 



In addition, by using (|3.6|) and the pre-smoothing property (|3.4p . respectively, 
for all x e C N ^ we obtain 

||GGG Fo,prca:||i < 0\\V O 

\\V0, P reX\\ 2 AoD -i Ao < a-\(\\x\\ 2 Ao - \\V , pTe x\\ Ao ). 



Hence 



Jjp\\CGC V , pTe x\\ 2 Ao < |M|^-||VWH^ 

< \\x\\ 2 Ao - \\CGC Vo, pTe x\\ 2 Ao , 
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since 

\\CGC Vo, pie x\\ Ao < \\CGC f Ao \\V , pie x\\ 2 Ao = \\V , pie x\\ 2 Ao , 
being \\CGC Q \\ An = 1. Therefore, for all x £ C N ^ n °\ it holds 

\\CGC V , pre xf Ao <(l + ^j 1 ||*||^. (A.4) 

By using inequality (IA.4|) in (|A.3[) , we have 

\\V 0tPOSt CGC V 0>pie x\\X < ]~y st/ ^ \\x\\ Ao , 

1 -+- CKprc/ p 

and the proof of Theorem 13.21 is concluded. 
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